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Abstract 

We present a new algorithm for solving an eigenvalue problem for a real sym- 
metric arrowhead matrix. The algorithm computes all eigenvalues and all com- 
ponents of the corresponding eigenvectors with high relative accuracy in 0{n^) 
operations. The algorithm is based on a shift-and-invert approach. Double 
precision is eventually needed to compute only one element of the inverse of 
the shifted matrix. Each eigenvalue and the corresponding eigenvector can be 
computed separately, which makes the algorithm adaptable for parallel com- 
puting. Our results extend to Hermitian arrowhead matrices, real symmetric 
diagonal-plus-rank-one matrices and singular value decomposition of real trian- 
gular arrowhead matrices. 
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1. Introduction and Preliminaries 

In this paper we consider eigenvalue problem for a real symmetric matrix 
A which is zero except for its main diagonal and one row and column. Since 
eigenvalues are invariant under similarity transformations, we can symmetrically 
permute the rows and the columns of the given matrix. Therefore, we assume 
without loss of generality that the matrix Ais anxn real symmetric arrowhead 
matrix of the form 
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where 



D = diag(di,(i2,...,d„-i) 
is diagonal matrix of order n — 1 , 



2 



(2) 



is a vector and a is a scalar. 

Such matrices arise in the description of radiationless transitions in isolated 
molecules Q , oscillators vibrationally coupled with a Fermi liquid [8| , quantum 
optics [iBl (see also Example S]). Such matrices also arise in solving symmetric 



real tridiagonal eigenvalue problems with the divide- and-conquer method 11 1. 

In this paper we present an algorithm which computes all eigenvalues and 
all components of the corresponding eigenvectors with high relative accuracy in 
O(n^) operations. 

Without loss of generality we may assume that A is irreducible, that is, 



and 



Ci^O, for aU 



di 7^ dj , for all i ^ j, i, j ~ I, . . . ,n — I. 



If A has a zero in the last column, say Ci = 0, then the diagonal element di 
is an eigenvalue whose corresponding eigenvector is the i-th unit vector, and 
we can reduce the size of the problem by deleting the i-th row and column of 
the matrix, eventually obtaining a matrix for which all elements Q are nonzero. 
If di — dj, then di is eigenvalue of matrix A (this follows from the interlacing 
property (O), and we can reduce the size of the problem by annihilating Q with 
a Givens rotation in the («,j)-plane and proceeding as in the previous case. 

Further, by symmetric row and column pivoting, we can order elements of 
D such that 

di > d2 > ■ ■ ■ > dn^i. (3) 

Hence, we will consider only ordered and irreducible arrowhead matrices. With- 
out loss of generality we can also assume that Ci > for *> which can be at- 
tained by pre- and post-multiplication of the matrix A with D = diag(sign(Ci))). 
Let 

A = VAV^ (4) 
be the eigenvalue decomposition of A. Here 

A = diag(Ai, A2, . . . , A„) 

is a diagonal matrix whose diagonal elements are the eigenvalues of A, and 
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is an orthonormal matrix whose columns are the corresponding ei genv ectors. 
The eigenvalues of A are the zeros of the Pick function (see [i, 16|) 



/(A)=a-A-^ 



i=l 



a- X- z^{D - Xiyh, 



(5) 



and the corresponding eigenvectors are given by 



Vi 



-1 



(6) 



Diagonal elements of the matrix D, di, are called poles of the function /. 

Notice that ([3]) and the Cauchy interlacing theorem [loj, Theorem 8.1.7] 
applied to matrices D and A imply the interlacing property 



Ai > di > A2 > ^2 > • • • > dn-2 > A„_i > dn-l > A^ 



(7) 



Since A is symmetric, its eigenvalues may be computed by invoking any 
of a number of standard programs (LAPACK However, these programs 

usually begin with an initial reduction of the matrix to tridiagonal form [17], or 
as proposed in [l6| . with an alternative which takes advantage of the structure 
of A by finding the zeros of the Pick function given in ([5]) , for the eigenvalues 
of A. This results in an algorithm which requires only 0{n^) computations 
and 0{n) storage. Although the idea is conceptually simple and in fact has 
been used to solve other eigenvalue problems of special structure ^, i5|, 1^, 0] , 
the computation is not always stable ll|. Namely, if the computed eigenvalues 
A,; are not accurate enough, then the computed eigenvectors Vi may not be 
sufficiently orthogonal (see Example [3]). The existing algorithms for arrowhead 
matrices (ill [l6j obtain orthogonal eigenvectors with the following procedure: 

- compute the eigenvalues Xi of A by solving ([5]); 



construct a new matrix 



A 



D 

~T 

z 



z 
a 



by solving inverse problem with the prescribed eigenvalues A, and diagonal 
matrix D, that is, compute new z and d as 



\ 



. w- ^ ■ (a, -ill) (Aj-ii,) 



J=2 
n-1 



= A„ + ^ [Xj - dj 



compute eigenvectors of A by 
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Since the formulas for Q involve only multiplications, division and subtrac- 
tions of exact quantities, each Q is computed with relative error of 0(£m)j 
where Sm denotes the machine precision^ Therefore, A — A + SA, where 
||(5^||2 = 0{eM)- Here |1 • II2 denotes the spectral matrix norm. We conclude 
that the computed eigenvalues satisfy standard perturbation bounds like 
those from [10|, Corollary 8.1.6]. Further, since \i are the eigenvalues of the 
matrix A computed to higher relative accuracy, the eigenvectors computed by 



([S]) are orthogonal to machine precision. For details see [11, 

Our algorithm uses a different approach. Accuracy of the eigenvectors and 
their orthogonality follows from high relative accuracy of the computed eigen- 
values and there is no need for follow-up orthogonalization. The algorithm is 
based on shift-and-invert technique. Basically, the eigenvalue A is computed as 
the largest or the smallest eigenvalue of the inverse of the matrix shifted to the 
pole di which is nearest to A, that is. 



A^i 

V 



(8) 



where v is either smallest or largest eigenvalue of the matrix 



a: 



[A -day 



Inverses of arrowhead matrices are structured in the following manner (here 
X stands for non-zero element): the inverse of an arrowhead matrix with zero 
on the shaft is a permuted arrowhead matrix with zero on the shaft, 



X 
X 

X 

X X 
XXX 



X 
X 

XXX 
X X 

X 



and the inverse of the full arrowhead matrix is a diagonal-plus-rank-one (DPRl) 
matrix, 



X X 

X X 

X X 

X X 

X X X X X 



± uu 



•^The machine precision ej\,/ is defined as a smallest positive number such that in the 
floating-point arithmetic 1 + em ¥^ 1- In Matlab or FORTRAN REAL(8) arithmetic cm = 
2.2204 ■ 10~^®, thus the floating-point numbers have approximately 16 significant decimal 
digits. The term "double of the working precision" means that the computations are performed 
with numbers having approximately 32 significant decimal digits, or with the machine precision 
equal to ej^^. 
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Our algorithm is completely parallel, since the computation of one eigen- 
value and its eigenvector is completely independent of the computation of other 
eigenvalues and eigenvectors. 

In Section 2 we describe the basic idea of our algorithm named aheig (Ar- 
rowHead EIGenvalues) . In Section 3 we discuss the accuracy of the algorithm. 
In Section 4 we present the complete algorithm which uses double of the working 
precision, if necessary. In Section 5 we illustrate algorithm with few examples 
and in Section 6 we apply our results to eigenvalue decomposition of Hermitian 
arrowhead matrix, singular value decomposition of real triangular arrowhead 
matrix and eigenvalue decomposition of real symmetric diagonal-plus-rank-one 
matrix. The proofs are given in Appendix A. 



2. Basic shift-and-invert algorithm 

Let A be an eigenvalue of A, let v be its eigenvector, and let x be the 
unnormalized version of v from Let di be the pole which is closest to A. 
Clearly, from Q it follows that cither A = Ai or A = Ai+i. Let Ai be the shifted 
matrix 

' Di 





where 



Ai = A- dj = 





D2 



Zl 
Z2 

a 



(9) 



Di = diag(di ~ di,. 
D2 = diag(di+i - di 

^1 = [ Ci C2 • • • 

Z2 = [ Q+1 Ci+2 

a = a — dj. 



■ , di-i — di), 

. . . , dn— 1 di) , 

Cn-l 1 , 



Notice that Di {D2) is positive (negative) definite. 
Obviously, if A is an eigenvalue of A, then 

II = X — di 

is an eigenvalue of Ai, and vice versa, and they both have the same eigenvector. 
The inverse of Ai is 



AJ 







Wi 








1 _ 


T 


b 


T 
W2 


1/G 







W2 


D2' 



















(10) 
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where 



W2 



Si 



Notice that 
where 



b = ^{-a + zfD^ ^zi + z^D^ ^22) 



f{d,)=a-d,-z^ [D-djy 



(11) 



where D is the diagonal matrix D without di and z is z without Q. 
The eigenvector x from ([6]) is given by 

[Di - fiiy^ Zl 

{D2 - III) Z2 
-1 



(12) 



If A is an eigenvalue of A which is closest to the pole di 
eigenvalue of matrix Ai which is closest to zero and 



then \i is the 



1 



± U; 



In this case, if all entries of A~ are computed with high relative accuracy, 
then, according to standard perturbation theory, v is computed to high relative 
accuracy (by any reasonable algorithm). In Section 3 we show that all entries 
of are indeed computed to high relative accuracy, except possibly h (see 
pT|) ). If h is not computed to high relative accuracy and it influences ||^i~^||2> 
it is sufficient to compute it in double of the working precision (see Section 4). 

Further, if /i is not the eigenvalue of Ai which is closest to zero, then \v\ < 
||^~^|L, and the quantity 



if. - 



(13) 



tells us how far is u from the absolutely largest eigenvalue of A^ . If K^, ^ 1 , 
then the standard perturbation theory does not guarantee that the eigenvalue 
II will be computed with high relative accuracy. Remedies of this situation are 
described in Remark |3l 

With this approach the componentwise high relative accuracy of the eigen- 
vectors computed by (fT2)) follows from high relative accuracy of the computed 
eigenvalues (see Theorem[3]). Componentwise high relative accuracy of the com- 
puted eigenvectors implies, in turn, their orthogonality. 
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The described procedure is implemented in algorithm aheigjbasic (Algo- 
rithm 1). The computation of the inverse of the shifted matrix, A~^, according 
to formulas PU)) and (|lip . is implemented in Algorithm 2. Algorithm 3 com- 
putes the largest or the smallest zero of the Pick function ^ by bisection. 
Given eigenvalue A, Algorithm 4 computes the corresponding eigenvector by ^ 
or respectively. 



3. Accuracy of the algorithm 



We now consider numerical properties of Algorithms 1, 2, 3, and 4. We 
assume tha standard model of floating_point arithmetic where subtraction is 



preformed with guard digit, such that 
fl(aob) = (ao6)(l + £o), 



e {+,-,*,/}, 



where Sm is machine precision. In the statements of the theorems and their 
proofs we shall use the standard first order approximations, that is, we neglect 
the terms of order 0{e\j) or higher. Moreover, we assume that neither overfiow 
or underflow occurs during the computation. 
We shall use the following notation: 



Matrix 


Exact eigenvalue 


Computed eigenvalue 


A 


A 


A 


A^ 


M 




A, = //(AO 




= IKV-) 


A,r^ 






{A-') = fliA-') 


v 





(14) 



Here 



A^ = fl (AO = 



Di {I + Ei) 














D2 (I + E2 



Zl 
Zl 

a(l + £a) 



where E\ and E2 are diagonal matrices whose elements are bounded by Em in 
absolute values and \ea\ < em- 

Further we define the quantities ka, and Kh as follows: 



h 

We also define the quantity 



-//(A) = A(H 


- kxEm) , 


= // (m) -A^(H 




= // (5) = /,(!-!- 


KfcEM) • 


\a\ + zfDj-^Zi 


+ \z^D-^Z2\ 


|— a + zfD^^zi 


+ zTD^^Z2\ 



(15) 
(16) 
(17) 



(18) 
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Algorithm 1 

[A, v] = aheig_basic {D, z, a, k) 

% Computes the fc-th eigenpair of an irreducible arrowhead matrix 
% A = [diag {D) z] z' a\ 
n = max{size{D)) + 1 

% Determine the shift tr, the shift index i, and whether A is on the left 
% or the right side of the nearest pole. 
% Exterior eigenvalues (fc = 1 or = n): 
if fc == 1 

a = di 

i = 1 

side = 'R' 
elseif k == n 

a = dn-i 

i = n — I 

side = 'L' 
else 

% Interior eigenvalues {k € {2, . . . , n — 1}): 

Dtemp — D — dk 
atemp = a — dk 
middle = Dte'mpk-i/2 

Fmiddle = atemp — middle — '^{z^-/ {Dtemp — middle)) 
if Fmiddle < 

a = dk 

i = k 

side = 'R' 
else 

a = dk-i 

i = k-l 

side = 'L' 
end 
end 

% Compute the inverse of the shifted matrix, 

[invDi,invD2, Wi,W2,Wt^,b] = invA(I?, z, a, i) 

% Compute the leftmost or the rightmost eigciivahic of A^^ 

V = bisect ([mt;Di; 0; invD2\, [wi; w^; W2\,b, side) 
% Compute the corresponding eigenvector 

jj. = l/v 

V = vect{D — (7, z, fi) 

% Shift the eigenvalue back 
X = fi + a 
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Algorithm 2 



[invDi,invD2T'Wi,W2,'W(^, b] — invA {D, z, a, i) 

% Computes the inverse of an arrowhead matrix A = [diag(D — di) z; z' a — di] 
% according to JlO} and (fTTj) . 

n ~ max(size(D)) + 1 
D = D-d, 
a ~ a — di 

wi = -zi.,i^i./Di.,i^i/zi 

W2 = —Zi+l:n-l./ Di^i;n-l/ Zi 

invDi = l./^Diyi-i 
invD2 = 

b = (-a + SMTO(2:i:j_l."2./Dl:i_l) + SMm(2;i+i:„_i.~2./Di+i:„_i))/zr2 



Algorithm 3 

A = bisect {D, z, a, side) 

% Computes the leftmost (for sjde='L') or the rightmost (for side—^K') eigenvalue 
% of an arrowhead matrix A = [diag (D) z; z' a] by bisection. 
n — max{size{D)) + 1 

% Determine the starting interval for bisection, [left, right] 
if side =— 'L' 

left = min{_D — \z\,a — \\z\\i} 

right = min di 
else 

right = maxjZ) + |z|,q:+ ||2;||i} 
left — laaxdi 
end 

% Bisection 

middle — (left + right) /2 
while {right — left)/abs{middle) > 2 * eps 
Fmiddle — a — middle — sum{z.~2./ (D — middle)) 
if Fmiddle > 
left — middle 
else 

right = middle 
end 

middle — {left + right) /2 
end 

% Eigenvalue 
A = right 
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Algorithm 4 

V = vect {D, z, A) 

% Computes the eigenvector of an arrowhead matrix A = [diag(£)) z; z' a] 
% which corresponds to the eigenvalue A by using (|6]). 
v=[z./{D-\):^l] 

V = v/\\v\\2 



3.1. Connection between accuracy of X and ^ 
Let 

X = fi + di 

be an eigenvalue of the matrix A, where fj, is the corresponding eigenvahie of 
the shifted matrix A^ = A — di from which A is computed. Let 

X = fl{Jl + d,) 

be the computed eigenvalue. Theorem [1] gives us dependency of accuracy of A 
in ([T5l) upon accuracy of Jl in ([T6| . 



Theorem 1. For A and X from I115\) and and /i from f we have 

|«a|<^^^(KI + 1). (19) 

Proofs of this theorem and subsequent theorems are given in Appendix A. 
From Theorem [1] we sec that the accuracy of A depends on and the size 
of the quotient 

M±M. (.0, 

Theorem [5] analyzes the quotient ([20)1 with respect to the position of A and signs 
of /i and the neighboring poles. 

Theorem 2. Let the assumptions of Theorem]^ hold, 
(i) If (see FigureUl (i)) 

sign(di) = sign(/i) , 

then 

\di\ + \^l\ ^ . 

|A| 

(ii) If A is between two poles of the same sign and s\gn{di) ^ sign(/i) (see 
Figure Ul (ii)), then 

\d^\ + \^^\ , „ 

Theorem [5] does not cover the foUowing cases: 
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di+i X di di X di-i 

(i) (ii) 

Figure 1: Typical situations from Theorem [2] 

(a) If di < 0, then ^ > 0. If, further, « \fJ,\, then Ai is near zero, and 
(|di| + ImI)/ I All » 1 (see Figure [l(a))g 

(b) If dn > 0, then /i < 0. If, further, \dn\ ~ then A„ is near zero, and 
again (|d„| + |Ai|)/|A„|» 10 

(c) If A is between two poles of the different signs and sign (di) ^ sign (/i), then 
either dij^\ < < di and < 0, or c?i < < o?i_i and fi > 0. In both cases, 
if, additionahy, \di\ « then A is near zero, and {\di\ + |/i|)/ |A| ^ 1 (see 
Figured (c)). 



\ • \ \ \ • \- 

di Ai di+i OA di 

(a) (c) 

Figure 2: Typical situations for special cases 

Since only one of these three cases can occur. Theorems [T] and [5] imply that 
for all eigenvalues X & a {A) , but eventually one, it holds 

M.I + ImI < 3^ 



|A| 

If one of the above cases does occur, remedies are given in the following 
remark. 

Remark 1. If one of the cases (a), (b) or (c) occurs, then A is an eigenvalue 
of A nearest to zero, and we can accurately compute it from the inverse of A. 
Notice that the inverse is of an unreduced arrowhead matrix with non-zero shaft 



In this case Ai is computed as a difference of two close quantities and cancellation can 
occur. 

^In this case is computed as a difference of two close quantities and cancellation can 
occur. 
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is a diagonal-plus-rank-one (DPRl) matrix of the form 



A- 







puu 



where 



Eigenvalues of A ^ are zeros of (see 

n 

^(A) = l + p^ 



a-zTD-^z 



j=i J 



X' 



Since the absolutely largest eigenvalue of is computed accurately according 
to standard perturbation theory, and 1/|A| = ||yl~^||2, A is also computed with 
high relative accuracy. In computing matrix , eventually p needs to be com- 
puted in higher precision. For more details see Remark [31 If the denominator 
in p is computed as zero, the matrix A is numerically singular and we can set 
A = 0. Notice that all components of the the corresponding eigenvector are still 
computed accurately. 

Remark 2. Notice that Algorithm 1 (and, consequently. Algorithm 5 below) 
can be easily modified to return both quantities, di and p, such that X = di+ p. 
If none of the remedies from Remark [1] were needed, these two quantities give 
additional information about A (that is, they give a more accurate representation 
of A) . An example is given in Example [2] 

We still need to bound the quantity from ()19|) . This quantity essentially 
depends on the accuracy of fl{b). The bound for is given in Theorem |6l 

3.2. Accuracy of the eigenvectors 

Since the eigenvector is computed by (IT2]) . its accuracy depends on the ac- 
curacy of p as described by the following theorem: 

Theorem 3. Let U6\) hold and let 



Xi 



= fK 



{Di {I + Ei)-piy'zi 

_g 

-1 



(21) 



be the computed un-normalized eigenvector corresponding to p and X. Then 
Xj ^ Xj {l + Ex^) , |£xj < 3(|k,J -|-3)eM, j = l,...,n. 

In other words, if is small, then all components of the eigenvector are 
computed to high relative accuracy. Since the accuracy of A and x depends 
on the accuracy of p (on the size of k^) in the next three subsections tells we 
discuss the accuracy of p. Since p is computed as an inverse of the eigenvalue 
of the matrix fl{A^^), we first discuss the accuracy of that matrix. 
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3.3. Accuracy of the matrix A^ ^ 
We have the following theorem: 

Theorem 4. For the computed elements of the matrix A^^ from U0\) and II 11]) 
for all (j, k) 7^ we have 

For the computed element b = (A^^)^. from |_?9| ) we /laue 

l^bl < (n + 3)i^b, 

where Kf, is defined by il8\) . 



The above theorem states that all elements of the matrix A^ are computed 
with high relative accuracy except possibly b. Therefore, we have to monitor 
whether b is computed accurately, and, if not, it needs to be computed in double 
of the working precision (see Section 4 for details). 

3.4. Accuracy of bisection 

Let Amax be the absolutely largest eigenvalue of a symmetric arrowhead 
matrix A, an let Amax be the eigenvalue computed by bisection as implemented 
in Algorithm 3. The error bound from 

0, 

Section 3.1] immediately implies 

that 

Amax Amax 

|T \ = KbisSM, Hbis < l-06n (%/n + l) . (22) 

I ^max I 

Notice that the similar error bound holds for all eigenvalues which are of the 
same order of magnitude as |Amax|- 

3.5. Accuracy of exterior eigenvalues of A^^ 

The desired interior eigenvalue and, in some cases, also absolutely smaller 
exterior eigenvalue A of A is in Algorithm 1 computed by where ly is one of 
the exterior eigenvalues of the matrix A~^. 

The following theorem covers the case when ly is the absolutely largest eigen- 
value of II A"^!!^, and gives two different bounds. 

Theorem 5. Let A^^ be defined by U0\) and let v be its eigenvalue such that 

\v\^\\A-^\\^. (23) 

Let V be the exact eigenvalue of the computed matrix i^A'^^^ — fl (A,^^). Let 

+ k^em) ■ (24) 

Then 



( 2 ""^ 

\k,\ < min (n + 3)^^i^6, 3^^ + (n + 3) (l + — ^ |Cfe|) 

I lvi| 4.-1 

where Kf, is defined by il8\) . 



k=l 



(25) 
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3.6. Final error hounds 

All previous error bounds are summarized as follows. 

Theorem 6. Let A he the computed eigenvalue of an unreduced arrowhead ma- 
trix A, let Jl he computed eigenvalue of the matrix Ai from (0), and let v he the 

corresponding computed eigenvalue of the matrix (^4^^^) from If ^ is the 

eigenvalue of Ai closest to zero ( or, equivalently, if i23\) holds ), then the error 
in the computed eigenvalue X is given hy U5\) with 

|ka| <3(|K,| + Kf„,)+4, (26) 

and the error in the computed un-normalized eigenvector x is given hy Theorem 
\^with 

\nA<K\+Kus + l. (27) 
where is bounded by ^2S\) and K^is is defined hy I12S\) . 

Since we are essentially using the shift-and-invert technique, we can guaran- 
tee high relative accuracy of the computed eigenvalue and high componentwise 
relative accuracy of the computed eigenvector if v is such that \v\ = 0{\\A~^\\2) 
and it is computed accurately. This is certainly fulfilled if the following condi- 
tions are met: 

CI. The quantity from \1S^) is moderate, and 

C2. (i) either the quantity Ki, from il8\) is small, or 

(ii) the quantity tt-t ^ from !i25\) is of order 0{n). 

k=l 

The condition CI implies that v will be computed accurately according to 
the standard perturbation theory. The conditions C2 (i) or C2 (ii) imply that 
Ki, from (j25l) is small, which, together with CI, implies that v is computed 
accurately. 

If the condition CI does not hold, that is, if ^ 1, remedies are given 
in Remark 2 below. If neither of the conditions C2 (i) and C2 (ii) holds, the 
remedy is to compute b in double of the working precision as described in Section 
4. 

Remark 3. We have two possibilities: 

(a) we can compute A by shifting to another neighboring pole provided that 
Ku is in this case small (shifting to the pole instead of di in Figure 
131(a)), 

(b) if shifting to another neighboring pole is not possible {Ki, ^ 1, see Figure 
13] (b)), we can invert A — al, where shift a is chosen near A, and a ^ 
{A, di, di_i}. This results in a DPRl matrix 



{A-aiy 



{D - aiy 



T 

pun , 
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where 



a-zT{D-aI)-^z' 
Eigenvalues of this matrix are zeros of 



n 2 



and the absolutely largest eigenvalue is coniputed accurately. Eventually, 
p needs to be coniputed in higher precision^ 



di di-1 di 

• \ — ^ • h- 

K+i A(Ai) Ai_i Aj+i A(Ai) Ai_i 

(a) (b) 

Figure 3: Typical situations from Remark [3] 



4. Final algorithm 

If neither of the conditions C2 (i) and C2 (ii) hold, in order to guarantee 
that A will be computed with high relative accuracy, the element h from the 
matrix needs to be computed in higher precision. The following theorem 
implies that if 1 <C Kb < 0{l/eM), it is sufficient to evaluate ([TT]) in double of 
the working prccision01f| 

Theorem 7. // -a > in if77]) . set 

P = -a + z'^D^^zi, Q = -z'^D^^Z2, 
and if —a < in ill]) set 

P = z^D^'^zi, Q ^a- z^D^^Z2. 



"Determining whether p needs to be computed in higher precision is done similarly as 
determining whether element b of needs to be computed in higher precision, which is 

described in Section 4. Further, Theorem [7] implies that it suffices to compute p in double of 
the working precision. 

''^If K), > 0{1/em), that is, if K), = l/es ior some se < £M, then, in view of Theorem [7] 
b needs to be computed with extended precision £e- 

* Usage of higher precision in conjunction with the eigenvalue computation for DPRl ma- 
trices is analyzed in jj], but there the higher precision computation is potentially needed 
in the iterative part. This is less convenient than our approach where the higher precision 
computation is used only to compute one element. 
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Notice that in both cases P,Q > and b = (P - Q)/C?- Let P = fl{P) and 
Q = fKQ) ^6 evaluated in standard precision, sm- Assume that P ^ Q and 
Kb < 0{\/eM)- If P, Q and b are all evaluated in double of the working 
precision, e\,j, then |_?7| ) holds with < 0{n). 

We summarize the above results in one, complete algorithm, aheig. The 
algorithm first checks the components of the vector z. If they are of the same 
order of magnitude, the eigenpair (A, v) is computed by Algorithm 1. If that is 
not the case, the quantity Ki, is computed, and if 3> 1, the eigenpair (A, u) 
is computed by Algorithm 1 but with evaluation of b in double of the working 
precision. At the end, the quantity K^, is computed, and if ^ 1, one of the 
remedies from Remark |3] is applied. 



Algorithm 5 

[A, v] = aheig {D, z, a, k) 

% Computes the k-th eigenpair of an ordered irreducible arrowhead matrix 
% A = [diag {D) z; z' a] 

compute the shift i as in the first part of Algorithm 1 
if the quantity (^J2 lOI^/ICil from ((25} is of 0{n) 

% standard precision is enough 
[A, v] = aheig_basic(i?, z, a, k) 
else 

compute the quantity Kjj from (1181) 
if i^f, > 1 
% double precision is necessary 

[A, v] ~ aheig_basic(_D, z, a, k) with evaluation of b in double precision 
else 

% standard precision is enough 
[A, v] = aheig_basic(D, z, a, k) 
end 
end 

compute the quantity from (I13p 
if i^^ > 1 

apply one of the remedies from Remark [3] 
end 



4-1. On implementing double precision 

Implementation of the double of the working precision depends upon whether 
the input is considered to be binary or decimal. 

Double standard precision in Matlab, which assumes that input is binary, 
is obtained by using a combination of commands vpa, digits and double |l3|, 
where 
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- digits (d) specifies the number of significant decimal digits d used to do 
variable precision arithmetic vpa, 

- vpa(x) uses variable-precision arithmetic to compute x to d decimal digits 
of accuracy, 

- double (x) converts x to standard precision. 

The assignment al=vpa(a,32) pads the binary representation of a with 
zeros, which means that the decimal interpretation of the variable al may have 
non-zero entries after 16-th significant decimal digit. The same effect is obtained 
in Intel FORTRAN compiler if ort [l^l by the following program segment 

real(8) a 
real(16) al 

al=a 

However, the user can assume that the true input is given as a decimal 
number, which is, for example, assumed by extended precision computation in 
Mathematica 20] . In this case, the options in Matlab are to either use symbolic 
computation, or to cast the input to a string, and then convert it to extended 
precision: 

al=vpa(num2str(a, 16) ,32) 

In this case, the the decimal interpretation of the variable al has all zero 
entries after 16-th significant decimal digit, but the binary representation of 
the variable a is, in general, padded with non-zero entries. The same effect 
is obtained in if ort writing to and reading from a string variable as in the 
following program segment: 

real(8) a 
real(16) al 
character (25) string 

writeCstring,*) a 
readCstring,*) al 

If the input consists of numbers for which decimal and binary representation 
are equal (for example, integers, as in Example [3] below) , then the two above 
approaches give the same results. 

5. Numerical Examples 

We illustrate out algorithm with four numerically demanding examples. Ex- 
amples [1] and [2] illustrate Algorithm 1 , Example [3] illustrates the use of double 
precision arithmetic, and Example |4] illustrates an application of higher dimen- 
sion. 
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Example 1. In this example both quantities K,y from and Kf, from (ITS)) are 
for all eigenvalues approximately equal to 1 , so we guarantee that all eigenvalues 
and all components of their corresponding eigenvectors are computed with high 
relative accuracy by Algorithm 5 (aheig) using only standard machine precision. 
Let 



A: 



2 • 10-3 














107 





10-7 











107 

















1 











-10-7 





107 














-2- 10-3 


107 


10^ 


10^ 


1 


107 


107 


1020 



The eigenvalues computed by Matlab [13| routine eig, Algorithm 5 and Math- 
ematica [l^ with 100 digits precision, are, respectively: 



1.000000000000000 • 10 
1.999001249000113 • 10 
4.987562099695390 • 10 
-1.000644853973479 • 10" 
-2.004985562101759 ■ 10 
-2.001001251000111 ■ 10 



20 
-3 
-9 
20 

-6 
-3 



1.000000000000000 • 10 
1.999001249000113 ■ 10 
4.987562099722817 ■ 10 
-9.999999999980001 ■ 10^ 
-2.004985562101717 ■ 10 
-2.001001251000111 ■ 10 



20 
-3 
-9 
20 

-6 
-3 



^(Math) 

1.000000000000000 ■ 10 
1.999001249000113 ■ 10 
4.987562099722817 • 10 
-9.999999999980001 ■ 10" 
-2.004985562101717 ■ 10 
-2.001001251000111 ■ 10 



20 
■3 
9 

20 

-6 
-3 



We see that even the tiniest eigenvalues A3 and A4, computed by Algorithm 5, are 
exact to the machine precision, which is not true for the eigenvalues computed 
by eig. Because of the accuracy of the computed eigenvalues, the eigenvectors 
computed by Algorithm 5 are componentwise accurate up to machine precision, 
and therefore, orthogonal up to machine precision. For example: 



,(e»9) 



(aheig) 



V 



(Math) 



4.999993626151683 ■ 


10" 


ll 


-4.999999999985000 


10" 


ll 


-4.999999999985000 


10" 


11 


9.999999962328609 


10" 


-7 


-9.999999999969000 


• 10 


-7 


-9.999999999969000 


■ 10 


-7 


9.999999999990000 


10" 


-1 


-9.999999999989999 


■ 10 


-1 


-9.999999999989999 


• 10 


-1 


-9.999999964673912 


• 10- 


-7 


9.999999999970999 


10" 


-7 


9.999999999970999 


10" 


-7 


-5.000006338012225 


10" 


11 


4.999999999985000 ■ 


10" 


11 


4.999999999985000 • 


10" 


11 


-9.999999963825105 


lo- 


21 


9.999999999970000 ■ 


lo- 


21 


9.999999999969999 


10"" 
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Example 2. In this example, despite very close diagonal elements, we again 
guarantee that all eigenvalues and all components of their corresponding eigen- 
vectors are computed with high relative accuracy, without deflation. Let 



A 



1 + 4:6 M 











1 





1 + Ssm 








2 








1 + 2eM 





3 











1 + l£M 


4 


1 


2 


3 


4 






where £m = 2 • 2 = 2.2204 • 10 For this matrix the quantities Ki, and 
Kb are again of order one for all eigenvalues, so Algorithm 5 uses only standard 
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working precision. The eigenvalues computed by Matlab and Algorithm 5 are: 



^(eig) y(aheig) 

6.000000000000000 6.000000000000001 

1 + 4eM 1 + 4eA/ 

1 + Aem 1 + 3£m 

1 + Sem 1 + 2£m 

-5.000000000000000 -4.999999999999999 

The eigenvalues computed by Mathematica with 100 digits precision, properly 
rounded to 32 decimal digits ar^^: 

6.0000000000000002018587317500285 
1.0000000000000008727792604471857 
1.0000000000000006206061701073114 
1.0000000000000003571862771540971 
-4.9999999999999998317843902083010 

The eigenvalues computed by Matlab are accurate according to standard per- 
turbation theory, but they do not satisfy the interlacing property. Furthermore, 
the Matlab's eigenvectors corresponding to A2, A3 and A4 only span an accurate 
eigenspace, and are not individually accurate. On the other hand, the eigen- 
values computed by Algorithm 5 are exact (they coincide with the eigenvalues 
computed by Mathematica properly rounded to 16 decimal digits). Notice that 
despite of very close eigenvalues. Algorithm 5 works without deflation. Due to 
the accuracy of the computed eigenvalues, the eigenvectors computed by Al- 
gorithm 5 are componentwise accurate up to the machine precision, and are 
therefore orthogonal. 

If, as suggested in Remark [2J the algorithms are modified to return di and 
/X (both in standard precision), then for the eigenvalues A2, A3 and A4 the cor- 
responding pairs {di , ^) give representations of those eigenvalues to 32 decimal 
digits. That is, exact values di + ^ properly rounded to 32 decimal digits are 
equal to the corresponding eigenvalues computed by Mathematica as displayed 
above. 

Example 3. In this example we can guarantee all eigenvalues and eigenvec- 
tors, componentwise will be computed with high relative accuracy only if we 
use double of the working precision when computing b from (jlip in matrices 



^Since, as described in Section 4.1, Mathematica uses decimal representation of the input, 
in order to obtain accurate eigenvectors we need to define em in Mathematica with the output 
of Matlab's command vpa(eps), em = 2.2204460492503130808472633361816 ■ lO-^^*. 
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\ An ^,A. ^ and A. ^ Let 



A 



The quantities K^, and if;, are: 






















4 











1 








3 








1 











2 





1 














1 


1 




1 


1 


1 


1 





10 



9.999999090793056 ■ 
9.999996083428923 • 

1.000000117045544 
9.999998561319470 • 

7.941165469988994 



10" 
10" 

10° 
10" 

10° 



3.243243243540540 ■ 10^ 
3.636363636818182 ■ 10^ 
4.444444445000000 ■ 10^ 
5.217390439488477 ■ 10^ 
5.217390439488477 ■ 10^ 



It is clear, from the condition numbers, that the element b in each of the ma- 



trices Ar, 



I ^3 ' 



A^ ^ and A^ ^ 



needs to be computed in double of the working 



)io_4 































1 








-1 








1 











-2 





1 














-3 


1 




1 


1 


1 


1 


IQio - 4 



precision. For example. 



A-d.I 



The element b = [A2^^]22 computed by Algorithm 2 gives b = 6.16666666667, 
Matlab's routine inv yields b = 6.166665889418350, while computing b in double 
of the working precision gives the correct value b = 6.166666668266667. 

Eigenvalues computed by Algorithm 1 {aheig_basic, using only standard 
working precision). Algorithm 5 {aheig, using double of the working precision to 
compute respective 5's) and Mathcmatica with 100 digits precision, respectively, 
are: 



\aheig.basic 








)^Math 




2.000000000000000 


lOio 


2.000000000000000 ■ 


101° 


2.000000000000000 ■ 


101° 


4.150396802313551 


■10° 


4.150396802279712 


■10° 


4.150396802279713 


■10° 


3.161498641452035 


•10° 


3.161498641430967 


■10° 


3.161498641430967 


■10° 


2.188045596352105 


■10° 


2.188045596339914 


■10° 


2.188045596339914 


■10° 


1.216093560005649 


■10° 


1.216093584948579 


■10° 


1.216093584948579 


■10° 


-7.160348702977373 • 


10-^ 


-7.160346250991725 ■ 


10-1 


-7.160346250991725 ■ 


10-1 



The eigenvectors computed by Algorithm 5 are componentwise accurate to ma- 
chine precision and therefore orthogonal. 



^"Algorithm 5 does not compute Ki, and Ki, for the first eigenvalue, since it is an absolutely 
largest one. 
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Example 4. This example comes from the research related to decay of excited 
states of quantmxi dots in in real photon crystals fl5|. In this case 

- a is quantum dot transition frequency, 

- di is a frequency of the z-th optical mode, and 

- is an interaction constant of the quantum dot with the i-th optical 
mode. 

The size of the matrix is changeable but, in realistic cases, it is between 10'^ and 
lO**. We ran a test example for n — 2501 where, typically, 

di e [5.87- 10", 1.38- 10^^], 

G e [1.05- 10^ 1.10- 10^], 

a = 9.7949881500060375-10". 

For this matrix the condition number Ki, ~ 1 for all eigenvalues and the compo- 
nents of the vector z do not differ by much in size, thus the conditions CI and 
C2 (ii) from Section 3 are fulfilled. Therefore, all eigenvalues and all compo- 
nents of all eigenvectors are computed with high relative accuracy by Algorithm 
5 using only standard working precision. On the other hand about half of the 
eigenvalues computed by the Matlab routine eig do not satisfy the interlacing 
property. 



6. Applications 

In this section we extend our results to eigenvalue decompositions of Her- 
mitian arrowhead matrices, singular value decompositions of real triangular ar- 
rowhead matrices and eigenvalue decompositions of real symmetric diagonal- 
plus-rank-one matrices. 



6.1. Hermitian arrowhead matrices 
Let 



D z 

z* a 



where 

D = diag((ii,d2, . . . ,(i„_i), 
is a real diagonal matrix of order n — 1, 

Z = [ Cl C2 ■ ■ ■ Cn-l ]* , 

is a complex valued vector and a is a real scalar. Here z* denotes the conjugate 
transpose of 2. As in Section 1, we assume that C is irreducible. The eigenvalue 
decomposition of C is given by 



C = UAU* 
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where A = diag(Ai, . . . , A„) S R"^" is a diagonal matrix of eigenvalues, and 
U = \ui U2 ■ ■ ■ Un] is an unitary matrix of the corresponding eigenvectors. 

To apply Algorithm 5 to Hermitian arrowhead matrix we first transform C 
to real symmetric arrowhead matrix A by diagonal unitary similarity: 



A = $*C$ 



where 



$ — diag 



ICil 



IC2I 



D \z\ 



Cn-l 



(28) 



■ ' ICn- 



1 



We now compute the k-th eigenpair (A, v) of A by Algorithm 5, and set u = 
$u. Since we guarantee high relative accuracy of the eigenvalue decomposition 
of A computed by Algorithm 5, we also guarantee high relative accuracy of the 
eigenvalue decomposition of C. Notice that, if double precision is needed to 
compute b in Algorithm 5, the modules |Ci| in (^5]) need to be computed in 
double of the working precision, as well. 

Remark 4. Similarly, for irreducible non-symmetric arrowhead matrix 



G 



D 

z 



z 
a 



where sign(^i) = sign(^j), i = l,...,n— 1, we define the diagonal matrix 



* = diag I sign(Ci)W ^, . . . , sign(C„-i) W y^, 1 
V Ci V C«-i 



The matrix 



A = -^-^G^ = 



D z 

T 

z a 



where G = V CiCi is an irreducible symmetric arrowhead matrix. 

We now compute the fc-th eigenpair (A, v) ol A by Algorithm 5. The eigenpair 
of G is then (A, set u — ^v. Since we guarantee high relative accuracy of 
the eigenvalue decomposition of A^ we also guarantee high relative accuracy of 
the eigenvalue decomposition of G. Notice that, if double precision is needed 
to compute h in Algorithm 5, the elements C,i need to be computed in double of 
the working precision, as well. 

6.2. Singular value decomposition of a triangular arrowhead matrix 
Let 



B = 



D 




z 
a 
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be an irreducible upper triangular arrowhead matrix, that is, di 7^ dj for i ^ j 
and Ci 7^ for all i. The matrix 



A = B^ B = 



D 



Dz 

a + z 



is an irreducible symmetric arrowhead matrix. 

When applying Algorithm 5 to the matrix A, we must ensure that all com- 
ponents of in (fTI?|) are computed to high relative accuracy. This is obviously 
true for elements of the vectors Wi and w^- Diagonal elements, except &, are 
computed with high relative accuracy as differences of squares of original quan- 
tities. 



The element h = [A';\ 



{d,-di){dj+di)' 
from (fTTj) is computed aJ"! 



6 = 



z + d,+ 



{dj - d,){dj + di 



If double precision is needed in Algorithm 5, all entries of A need to be computed 
in double precision. 

Let B = UYjV^ be the singular value decomposition of B, where E = 
diag(CTi , . . . , CT„) are the singular values, the columns of V are the corresponding 
right singular vectors and the columns of U are the corresponding left singular 
vectors. We first compute the fe-th eigenpair (A, v) of A by Algorithm 5. Then 
a = \/A is the corresponding singular value of B and v is the corresponding 
right singular vector. The value a and all components of v are computed to 
almost full accuracy. From the relation B = YjV'^ for the A:-th row we have 



T 
K.n- 



D z 
a 



which implies 



■Ul:„_l = (7Wl:„_lD ^ 

From the relation BV = C/S for the fc-th column we have 



which implies 



'D z 










a 






~ a 






Vn 




Un 



av„ 
(J 



Components of u are computed by multiplication and division of quantities 
which are accurate to almost full machine precision, so the are accurate to 
almost full machine precision, as well. 



^^In view of Theorem [T] if double precision computation is necessary, the positive and 
negative parts of this formula should be computed separately, and then added. 
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6.3. Diagonal-plus-rank- one matrices 
Let 



where 



D = diag(di, 
u = \ui 



M^D 



, dn), di> d2> ■ ■ ■ > dn, 
iT 



Ui 7^ 0, 1 = 1, 



be a n X n irreducible ordered real symmetric diagonal-plus-rank-one (DPRl) 
matrix. Let 

D = diag(di, . . . ,d„_i), 



U= [Ui ■ ■■ Un~l\ 

u„A-i O" 
-u^A-i 1 ■ 



L = 



Then 



where 



A = L-^ML = 



a 



D z 

T 

z a 



dn -\- U U. 



is an irreducible real symmetric arrowhead matrix. 

When applying Algorithm 5 to the matrix A, we must ensure that all com- 
ponents of in (fTO|) are computed to high relative accuracy. This is obviously 
true for elements of the vectors Wi and ^2- Diagonal elements, except 6, are 
computed with high relative accuracy as differences of original quantities, and 
the element b = from (jlip is computed as 



6 = T- I -dn - u^u + rf, + ^ 



If double precision is needed in Algorithm 5, all entries of A need to be computed 
in double precision. 

Let M = QAQ^ and A = VAV'^ be the eigenvalue decompositions of M 
and A, respectively. Since M is by assumption irreducible, its eigenvalues satisfy 
interlacing property 



Ai > di > A2 > ^2 > • • • > A„ > d„. 



(29) 



We first compute the fc-th eigenpair (A, v) of A by Algorithm 5. The value A 
and all components of v are computed to almost full accuracy. The relation 
V'^AV = V'^L~'^MLV = A implies that the columns of the matrix X = LV 
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are the unnormalized eigenvectors of the matrix M . Further, since, by (P^. ah 
eigenvalues are simple, we conclude that X = QE, where E — diag(cri, . . . , ct„) 
is a positive definite matrix. Notice that QYiV'^ — L is, in fact, singular value 
decomposition of L. 

Equating fc-th columns of the equation X ~ LV gives 





X 




■ u„A-i 0' 




V 


X — 










Vn_ 



where x and v are partitioned according to L. This immediately implies that 

X — UnA^^V. 

Notice that, since all components of v are computed to almost full, accuracy, the 
same holds for the components of x, and it remains to compute x„ accurately. 
Let 

' q 

be the fc-th column of Q and let a = E^fc. Equating fc-th rows of the equation 

X-^ = E-iQ^ = V^L-^ 

gives for the n-th element 



1 
a 



Thus, 



1 



and, in order to compute a;„ , it is necessary to compute . From X = UT, = LV 
it follows that V'^L'^LV = E^, or, equivalently, LV = L'^T/E^. Equating fc-th 
columns of this equation gives 



Av- 



1 



u v„ 



This gives n — 1 equations for ct^, and we can choose the numerically most 
accurate one. 

Therefore, Xn will be computed to almost full machine precision, as are the 
entries of x, and it remains to normalize x and obtain q = x/a. 

Remark 5. Notice that DPRl matrices of the form D — uvF cannot be reduced 
to symmetric arrowhead matrix by the procedure described in this section. By 
using ideas from this paper, it is possible to derive highly accurate algorithm 
for DPRl matrices without prior transformation to arrowhead form. This al- 
gorithm, which is a topic of our forthcoming paper, covers more general DPRl 
matrices of the form 



D 



puu 
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Appendix A. Proofs 

Proof of Theorem [TJ 

Let Jl and A be defined by (fHl) . Tlien 

X = fl {d, + = {di + A^) (1 + ei) • 
By simplifying the equality 

{di + fJ.{l + k^Em)) (1 + ei) = A (1 + kxEm) 
and using A = /i + d,; , we have 

diEi + H {Kf^EM + £i) = Aka£j\/- 
Taking absolute value gives 

M<^^(KH-i). □ 

Proof of Theorem 

(i) The assumption sign {di) ~ sign(/i) immediately implies 

\d^\ + \^^\ ^ \d^ + Ml ^ ^ 
|A| + 

(^iij The assumptions imply that either 

< di+i < X < di, fJ- <0, 

or 

di < A < < 0, fi> 0. 
In the first case A is closest to the pole di and 

\di\ + ^ Mil + ^ \di - di+i\ ^ di + ^di - \di+i 
|A| ~ \\di + di+i\ ~ + idj+i 

\di + jdi+i d'i 

Here we used the inequalities < \ \di — c?i+i| and |A| > 5 |di + d,;+i| for the 
first inequality, di — di+i > and d,; + dj+i > for the second inequality and 
di+i > for the fourth inequality, respectively. 
The proof for the second case is analogous. □ 
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Proof of Theorem [3 

Let X and x be defined by and (PTj) . respectively. The tlicorcm obviously 
holds for Xn = Xn = —i- For 5;^ we have 

x^ - f-i) = 7—^^ T (1 + £i) = X, (1 + . 

By using ([T5)) and (HH), the first order approximation gives 

k^J < (Ika^I + l)eA/. 
For J ^ by solving the equality 



" {{d, - d,) (1 + ei) - M (1 + A^pEAf)) (1 + £2) ^ " rf^A 



.(l + e3) = ^(l + e.) 

for £x, using (|16p and A = fi + di, and ignoring higher order terms, we have 
(dj - di) (ei + £2 + £3) - M (k/x£m + £2 + £3) 



£.T 



- A 



Therefore, 



|£.|<^^^^^^^(KI+3)£M. (A.1) 
To complete the proof we need to analyze two cases. If 
sign (dj - di) = - sign/^, 

then 

\dj - d^\ + 1^1 _ \dj - di- _ \dj - A| _ 
\d, - A| " Id, - A| ' \d, - Ar ■ 

If 

sign(dj - di) = sign^, 
then, since di is pole closest to A, we have < 0.5 \dj — di\ and 

\dj-d,\ + \^i\ ^ \dj-d,\ + \n\ ^ ^\d.j~di\ ^ ^ 
\dj-X\ - \d,-di\-\fi\- ^\dj-d,\ 

Finally, the theorem follows by inserting this into (|A.1|) . □ 
Proof of Theorem 

For the non-zero computed elements of the matrix A^^ from (fTO|) and (fTTjl . 
except the element b = [AY^]ii, we have: 

f^^^^'^'y-idTdWrir)^'^'^^' 

fli[Ar%.)^fl{[Ar%„)^hl + se), 
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where |£fe| < Em for all indices k. The first statement of the theorem now follows 
by using standard first order approximations. 
Similar analysis of the formula ([TT]) yields 

fl{[A-%) = b = b + Sb, 

where 

\Sb\ < ^ {\~a\ + \zfD^^zi\ + \ziD^^Z2\) {n + 3)eM- (A.2) 

This, in turn, implies (1171) with 

I 1 ^ Ja\+\zlD^'z,\ + \zjD^'z2\ 

\Kb\ < -rrr 1 r = + 3)-! tt^--\ — r = + 3)a6, 

\b\ \eM\ ' \-a + zfD^^zi+ z^D2^Z2\ 

where Kb is defined by (fTS]). □ 



Proof of Theorem O 

Let ^ 

Therefore, 

i9-H = iiMr^ii2, 

which, together with ([M)) . implies 

I^^k.emI < ll-JA-'lb- (A.3) 

Theorem m implies that 

\\SA-'\\2 < in + 3)\\\A;^\hKbeM- 
Since |!|Ari|||2 < y^\\Ar^\\2 and = ||Ari||2 , from JMl) we have 

K\ < {n + 3)V^Kb, (A.4) 



which proves the first part of the bound ([25)) . 

For the second part of the proof, notice that Theorem 2] also implies 

||Mri||2<3|||^|||2eM + |^fe|, (A.5) 

where A is equal to the matrix A^^ without 6 (that is, with An = 0). 

By bounding (|A.3|) with (|A.5P and (|A.2p . and dividing the resulting inequal- 
ity by I uEm I , we have 

< 3V«+ Ui + 3)[^|^ ■ ^2-^ ■ -J- (A.6) 
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Since 

\-a\ 1 



2 — /■2 I a + z'[D-^^Zi+ Z2D2^Z2- z'[D-^^Zi~ Z2D2^Z2\ 



Si s^ 

1 

c 



<\b\ + ^{\zfD~'z,\ + \z^D-'z2\ 



from (|A.6[) it follows 



< SV'T- + (n + 3) — + —- '—^ . (A. 7) 



Since |5| < \v\ and 



|A, ML li'l = max llA, ixll > llA^ ^efcl 
I i 112 ||^||2 = i" "2 II I 



a ^ \Ck\ 



> 



\{dk-d,f CHdk'd^f \Q\\dk~d, 
by simply dividing each term 



Cf\dk - d, 

in (jA.7| with the corresponding quotient 



ICil \dk - d^\ ' 

we obtain 



I'^.l < 3V^+(n + 3)(^l + -^jai j 



The bound 125]) now follows from (|A.4p and (jA.Sp . □ 



(A.8) 



Proof of Theorem O 

We first prove the bound (P7| . Since = /^(z?) is computed by bisection, 
from ([2H) we have 

This and imply 

? = + K^£a/)(1 + KUsSm)- 

Since Jl = fl(l/i)), the bound (|27l) follows by ignoring higher order terms. The 
bound (|26|) now follows by inserting ((27)) into Theorems [1] and [H □ 
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Proof of Theorem 

Let the assumptions of the theorem hold. Let b be computed in double of the 
working precision, e\j , and then stored in the standard precision. The standard 
floating-point error analysis with neglecting higher order terms gives 

-2 [i- + KiEj^^j) — — -2 (i + KfeeMj 

= 5(1 + Kb£A/) , 

where |kp|,|kq| < (n + 1) and |ki| < 3. Solving the above equality for k;,, 
neglecting higher order terms, and taking absolute values gives 

< ^^p~^qI + 4) £M = Kb{n + 4)£M. 

Since, by assumption, Ki, < 0(l/ej\/), this implies 

\Kb\ < 0(n), 

as desired. □ 
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